%%% Historical trends

clear
load('counter.mat')

%%% Figure 2:
cz_division = load('dsets_for_model\cz_division_cpp.csv');

decades = [1890:20:2010];

pop_growth_data = log(data_pop(:,3:end)) - log(data_pop(:,2:end-1));

lhs = pop_growth_data(:);
rhs = zeros(N*length(decades),length(decades)*9);
for t=1:length(decades)
    for i=1:N
        rhs((t-1)*N+i,(t-1)*9+cz_division(i,2)) = 1;
    end
end

betta_aux = rhs\lhs
lhs_res = lhs - rhs*betta_aux;
lhs_res_reshapes = reshape(lhs_res,N,length(decades));

figure
plot(decades,lhs_res_reshapes(159,:),'r','LineWidth',1,'MarkerSize',6,'MarkerEdgeColor','r','MarkerFaceColor','r','Marker','d')
ylabel('20-year growth rate of population')
hold on
plot(decades,lhs_res_reshapes(269,:),'b','LineWidth',1,'MarkerSize',6,'MarkerEdgeColor','b','MarkerFaceColor','b','Marker','s')
plot(decades,lhs_res_reshapes(414,:),'g','LineWidth',1,'MarkerSize',6,'MarkerEdgeColor','g','MarkerFaceColor','g','Marker','o')
plot(decades,zeros(1,length(decades)),'k--','HandleVisibility','off')

%%% Figure 1:
class8_crosswalk = [1 1 1 2 2 3 4 5 5 6 7];
n_classes = 7;
patent_shares = zeros(N,n_classes,T);
patent_shares_national = zeros(n_classes,T);

for t=1:T
   for i=1:N
       for s=1:S
           patent_shares(i,class8_crosswalk(s),t) = patent_shares(i,class8_crosswalk(s),t) + data_pat(i,s,t);
       end
       patent_shares(i,:,t) = patent_shares(i,:,t) / sum(patent_shares(i,:,t));
   end
   for s=1:S
       patent_shares_national(class8_crosswalk(s),t) = patent_shares_national(class8_crosswalk(s),t) + sum(data_pat(:,s,t));
   end
   patent_shares_national(:,t) = patent_shares_national(:,t) / sum(patent_shares_national(:,t));
end

decades = [1870:20:2010];
figure
subplot(2,2,1), area(decades,squeeze(patent_shares(159,:,2:end))')
ax = gca; % current axes
ax.XTick = decades;
title('Detroit','FontSize',14)

subplot(2,2,2), area(decades,squeeze(patent_shares(414,:,2:end))')
ax = gca; % current axes
ax.XTick = decades;
title('Austin','FontSize',14)

subplot(2,2,3), area(decades,squeeze(patent_shares(269,:,2:end))')
ax = gca; % current axes
ax.XTick = decades;
title('Boston','FontSize',14)

subplot(2,2,4), area(decades,squeeze(patent_shares_national(:,2:end))')
ax = gca; % current axes
ax.XTick = decades;
title('National','FontSize',14)

%%% Appendix Figure A.1:
class11_crosswalk = [1 2 3 4 5 6 7 8 9 10 11];
n_classes = 11;
patent_shares = zeros(N,n_classes,T);
patent_shares_national = zeros(n_classes,T);

for t=1:T
   for i=1:N
       for s=1:S
           patent_shares(i,class11_crosswalk(s),t) = patent_shares(i,class11_crosswalk(s),t) + data_pat(i,s,t);
       end
       patent_shares(i,:,t) = patent_shares(i,:,t) / sum(patent_shares(i,:,t));
   end
   for s=1:S
       patent_shares_national(class11_crosswalk(s),t) = patent_shares_national(class11_crosswalk(s),t) + sum(data_pat(:,s,t));
   end
   patent_shares_national(:,t) = patent_shares_national(:,t) / sum(patent_shares_national(:,t));
end

decades = [1870:20:2010];
figure
subplot(2,2,1), area(decades,squeeze(patent_shares(159,:,2:end))')
ax = gca; % current axes
ax.XTick = decades;
title('Detroit','FontSize',14)

subplot(2,2,2), area(decades,squeeze(patent_shares(414,:,2:end))')
ax = gca; % current axes
ax.XTick = decades;
title('Austin','FontSize',14)

subplot(2,2,3), area(decades,squeeze(patent_shares(269,:,2:end))')
ax = gca; % current axes
ax.XTick = decades;
title('Boston','FontSize',14)

subplot(2,2,4), area(decades,squeeze(patent_shares_national(:,2:end))')
ax = gca; % current axes
ax.XTick = decades;
title('National','FontSize',14)


XLabels = 1:S; YLabels = 1:S;
CustomXLabels = string(XLabels); CustomYLabels = string(YLabels);
CustomXLabels(1,1:11) = ["A1","A2","A3","B1","B2","C1","E1","F1","F2","G1","H1"];
CustomYLabels(1,1:11) = ["A1","A2","A3","B1","B2","C1","E1","F1","F2","G1","H1"];

%%% Appendix Figure A.6:
map = zeros(1000,3);
for i=1:1000
   map(1001-i,:) = [i/1000,i/1000,i/1000]; 
end
figure
data_aux = load('dsets_for_model\patents_data\citprob_classes_1950_1970.csv');
data_aux = reshape(data_aux,S,S)';
subplot(1,2,1), heatmap(data_aux, 'Colormap', map, 'XDisplayLabels', CustomXLabels, 'YDisplayLabels', CustomYLabels)
data_aux = load('dsets_for_model\patents_data\citprob_classes_1990_2010.csv');
data_aux = reshape(data_aux,S,S)';
subplot(1,2,2), heatmap(data_aux, 'Colormap', map, 'XDisplayLabels', CustomXLabels, 'YDisplayLabels', CustomYLabels)

%%% Appendix Figure A.5:
map = zeros(1000,3);
for i=1:1000
   map(1001-i,:) = [i/1000,i/1000,i/1000]; 
end
figure
data_aux = load('dsets_for_model\patents_data\citprob_classes_geq1950.csv');
data_aux = reshape(data_aux,S,S)';
heatmap(data_aux, 'Colormap', map, 'XDisplayLabels', CustomXLabels, 'YDisplayLabels', CustomYLabels)



